home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / randist / discrete.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-06-05  |  13.1 KB  |  395 lines
  1. /* randist/discrete.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /*
  21.    Random Discrete Events
  22.    
  23.    Given K discrete events with different probabilities P[k]
  24.    produce a value k consistent with its probability.
  25.  
  26.    This program is free software; you can redistribute it and/or
  27.    modify it under the terms of the GNU General Public License as
  28.    published by the Free Software Foundation; either version 2 of the
  29.    License, or (at your option) any later version.
  30.  
  31.    This program is distributed in the hope that it will be useful, but
  32.    WITHOUT ANY WARRANTY; without even the implied warranty of
  33.    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  34.    General Public License for more details.  You should have received
  35.    a copy of the GNU General Public License along with this program;
  36.    if not, write to the Free Foundation, Inc., 59 Temple Place, Suite
  37.    330, Boston, MA 02111-1307 USA
  38. */
  39.  
  40. /*
  41.  * Based on: Alastair J Walker, An efficient method for generating
  42.  * discrete random variables with general distributions, ACM Trans
  43.  * Math Soft 3, 253-256 (1977).  See also: D. E. Knuth, The Art of
  44.  * Computer Programming, Volume 2 (Seminumerical algorithms), 3rd
  45.  * edition, Addison-Wesley (1997), p120.
  46.  
  47.  * Walker's algorithm does some preprocessing, and provides two
  48.  * arrays: floating point F[k] and integer A[k].  A value k is chosen
  49.  * from 0..K-1 with equal likelihood, and then a uniform random number
  50.  * u is compared to F[k].  If it is less than F[k], then k is
  51.  * returned.  Otherwise, A[k] is returned.
  52.    
  53.  * Walker's original paper describes an O(K^2) algorithm for setting
  54.  * up the F and A arrays.  I found this disturbing since I wanted to
  55.  * use very large values of K.  I'm sure I'm not the first to realize
  56.  * this, but in fact the preprocessing can be done in O(K) steps.
  57.  
  58.  * A figure of merit for the preprocessing is the average value for
  59.  * the F[k]'s (that is, SUM_k F[k]/K); this corresponds to the
  60.  * probability that k is returned, instead of A[k], thereby saving a
  61.  * redirection.  Walker's O(K^2) preprocessing will generally improve
  62.  * that figure of merit, compared to my cheaper O(K) method; from some
  63.  * experiments with a perl script, I get values of around 0.6 for my
  64.  * method and just under 0.75 for Walker's.  Knuth has pointed out
  65.  * that finding _the_ optimum lookup tables, which maximize the
  66.  * average F[k], is a combinatorially difficult problem.  But any
  67.  * valid preprocessing will still provide O(1) time for the call to
  68.  * gsl_ran_discrete().  I find that if I artificially set F[k]=1 --
  69.  * ie, better than optimum! -- I get a speedup of maybe 20%, so that's
  70.  * the maximum I could expect from the most expensive preprocessing.
  71.  * Folding in the difference of 0.6 vs 0.75, I'd estimate that the
  72.  * speedup would be less than 10%.
  73.  
  74.  * I've not implemented it here, but one compromise is to sort the
  75.  * probabilities once, and then work from the two ends inward.  This
  76.  * requires O(K log K), still lots cheaper than O(K^2), and from my
  77.  * experiments with the perl script, the figure of merit is within
  78.  * about 0.01 for K up to 1000, and no sign of diverging (in fact,
  79.  * they seemed to be converging, but it's hard to say with just a
  80.  * handful of runs).
  81.  
  82.  * The O(K) algorithm goes through all the p_k's and decides if they
  83.  * are "smalls" or "bigs" according to whether they are less than or
  84.  * greater than the mean value 1/K.  The indices to the smalls and the
  85.  * bigs are put in separate stacks, and then we work through the
  86.  * stacks together.  For each small, we pair it up with the next big
  87.  * in the stack (Walker always wanted to pair up the smallest small
  88.  * with the biggest big).  The small "borrows" from the big just
  89.  * enough to bring the small up to mean.  This reduces the size of the
  90.  * big, so the (smaller) big is compared again to the mean, and if it
  91.  * is smaller, it gets "popped" from the big stack and "pushed" to the
  92.  * small stack.  Otherwise, it stays put.  Since every time we pop a
  93.  * small, we are able to deal with it right then and there, and we
  94.  * never have to pop more than K smalls, then the algorithm is O(K).
  95.  
  96.  * This implementation sets up two separate stacks, and allocates K
  97.  * elements between them.  Since neither stack ever grows, we do an
  98.  * extra O(K) pass through the data to determine how many smalls and
  99.  * bigs there are to begin with and allocate appropriately.  In all
  100.  * there are 2*K*sizeof(double) transient bytes of memory that are
  101.  * used than returned, and K*(sizeof(int)+sizeof(double)) bytes used
  102.  * in the lookup table.
  103.    
  104.  * Walker spoke of using two random numbers (an integer 0..K-1, and a
  105.  * floating point u in [0,1]), but Knuth points out that one can just
  106.  * use the integer and fractional parts of K*u where u is in [0,1].
  107.  * In fact, Knuth further notes that taking F'[k]=(k+F[k])/K, one can
  108.  * directly compare u to F'[k] without having to explicitly set
  109.  * u=K*u-int(K*u).
  110.  
  111.  * Usage:
  112.  
  113.  * Starting with an array of probabilities P, initialize and do
  114.  * preprocessing with a call to:
  115.  
  116.  *    gsl_rng *r;
  117.  *    gsl_ran_discrete_t *f;
  118.  *    f = gsl_ran_discrete_preproc(K,P);
  119.    
  120.  * Then, whenever a random index 0..K-1 is desired, use
  121.  
  122.  *    k = gsl_ran_discrete(r,f);
  123.  
  124.  * Note that several different randevent struct's can be
  125.  * simultaneously active.
  126.  
  127.  * Aside: A very clever alternative approach is described in
  128.  * Abramowitz and Stegun, p 950, citing: Marsaglia, Random variables
  129.  * and computers, Proc Third Prague Conference in Probability Theory,
  130.  * 1962.  A more accesible reference is: G. Marsaglia, Generating
  131.  * discrete random numbers in a computer, Comm ACM 6, 37-38 (1963).
  132.  * If anybody is interested, I (jt) have also coded up this version as
  133.  * part of another software package.  However, I've done some
  134.  * comparisons, and the Walker method is both faster and more stingy
  135.  * with memory.  So, in the end I decided not to include it with the
  136.  * GSL package.
  137.    
  138.  * Written 26 Jan 1999, James Theiler, jt@lanl.gov
  139.  * Adapted to GSL, 30 Jan 1999, jt
  140.  
  141.  */
  142.  
  143. #include <config.h>
  144. #include <stdio.h>              /* used for NULL, also fprintf(stderr,...) */
  145. #include <stdlib.h>             /* used for malloc's */
  146. #include <math.h>
  147. #include <gsl/gsl_rng.h>
  148. #include <gsl/gsl_randist.h>
  149. #define DEBUG 0
  150. #define KNUTH_CONVENTION 1      /* Saves a few steps of arithmetic
  151.                                  * in the call to gsl_ran_discrete()
  152.                                  */
  153.  
  154. /*** Begin Stack (this code is used just in this file) ***/
  155.  
  156. /* Stack code converted to use unsigned indices (i.e. s->i == 0 now
  157.    means an empty stack, instead of -1), for consistency and to give a
  158.    bigger allowable range. BJG */
  159.  
  160. typedef struct {
  161.     size_t N;                      /* max number of elts on stack */
  162.     size_t *v;                     /* array of values on the stack */
  163.     size_t i;                      /* index of top of stack */
  164. } gsl_stack_t;
  165.  
  166. static gsl_stack_t *
  167. new_stack(size_t N) {
  168.     gsl_stack_t *s;
  169.     s = (gsl_stack_t *)malloc(sizeof(gsl_stack_t));
  170.     s->N = N;
  171.     s->i = 0;                  /* indicates stack is empty */
  172.     s->v = (size_t *)malloc(sizeof(size_t)*N);
  173.     return s;
  174. }
  175.  
  176. static void
  177. push_stack(gsl_stack_t *s, size_t v)
  178. {
  179.     if ((s->i) >= (s->N)) {
  180.         fprintf(stderr,"Cannot push stack!\n");
  181.         abort();                /* FIXME: fatal!! */
  182.     }
  183.     (s->v)[s->i] = v;
  184.     s->i += 1;
  185. }
  186.  
  187. static size_t pop_stack(gsl_stack_t *s)
  188. {
  189.     if ((s->i) == 0) {
  190.         fprintf(stderr,"Cannot pop stack!\n");
  191.         abort();               /* FIXME: fatal!! */
  192.     }
  193.     s->i -= 1;
  194.     return ((s->v)[s->i]);
  195. }
  196.  
  197. static inline size_t size_stack(const gsl_stack_t *s)
  198. {
  199.     return s->i;
  200. }
  201.  
  202. static void free_stack(gsl_stack_t *s)
  203. {
  204.     free((char *)(s->v));
  205.     free((char *)s);
  206. }
  207.  
  208. /*** End Stack ***/
  209.  
  210.  
  211. /*** Begin Walker's Algorithm ***/
  212.  
  213. gsl_ran_discrete_t *
  214. gsl_ran_discrete_preproc(size_t Kevents, const double *ProbArray)
  215. {
  216.     size_t k,b,s;
  217.     gsl_ran_discrete_t *g;
  218.     size_t nBigs, nSmalls;
  219.     gsl_stack_t *Bigs;
  220.     gsl_stack_t *Smalls;
  221.     double *E;
  222.     double pTotal = 0.0, mean, d;
  223.     
  224.     if (Kevents < 1) {
  225.       /* Could probably treat Kevents=1 as a special case */
  226.  
  227.       GSL_ERROR_VAL ("number of events must be a positive integer", 
  228.             GSL_EINVAL, 0);
  229.     }
  230.  
  231.     /* Make sure elements of ProbArray[] are positive.
  232.      * Won't enforce that sum is unity; instead will just normalize
  233.      */
  234.  
  235.     for (k=0; k<Kevents; ++k) {
  236.         if (ProbArray[k] < 0) {
  237.       GSL_ERROR_VAL ("probabilities must be non-negative",
  238.                 GSL_EINVAL, 0) ;
  239.         }
  240.         pTotal += ProbArray[k];
  241.     }
  242.  
  243.     /* Begin setting up the main "object" (just a struct, no steroids) */
  244.     g = (gsl_ran_discrete_t *)malloc(sizeof(gsl_ran_discrete_t));
  245.     g->K = Kevents;
  246.     g->F = (double *)malloc(sizeof(double)*Kevents);
  247.     g->A = (size_t *)malloc(sizeof(size_t)*Kevents);
  248.  
  249.     E = (double *)malloc(sizeof(double)*Kevents);
  250.  
  251.     if (E==NULL) {
  252.       GSL_ERROR_VAL ("Cannot allocate memory for randevent", ENOMEM, 0);
  253.     }
  254.  
  255.     for (k=0; k<Kevents; ++k) {
  256.         E[k] = ProbArray[k]/pTotal;
  257.     }
  258.  
  259.     /* Now create the Bigs and the Smalls */
  260.     mean = 1.0/Kevents;
  261.     nSmalls=nBigs=0;
  262.     for (k=0; k<Kevents; ++k) {
  263.         if (E[k] < mean) ++nSmalls;
  264.         else             ++nBigs;
  265.     }
  266.     Bigs   = new_stack(nBigs);
  267.     Smalls = new_stack(nSmalls);
  268.     for (k=0; k<Kevents; ++k) {
  269.         if (E[k] < mean) {
  270.             push_stack(Smalls,k);
  271.         }
  272.         else {
  273.             push_stack(Bigs,k);
  274.         }
  275.     }
  276.     /* Now work through the smalls */
  277.     while (size_stack(Smalls) > 0) {
  278.         s = pop_stack(Smalls);
  279.         if (size_stack(Bigs) == 0) {
  280.             /* Then we are on our last value */
  281.             (g->A)[s]=s;
  282.             (g->F)[s]=1.0;
  283.             break;
  284.         }
  285.         b = pop_stack(Bigs);
  286.         (g->A)[s]=b;
  287.         (g->F)[s]=Kevents*E[s];
  288. #if DEBUG
  289.         fprintf(stderr,"s=%2d, A=%2d, F=%.4f\n",s,(g->A)[s],(g->F)[s]);
  290. #endif        
  291.         d = mean - E[s];
  292.         E[s] += d;              /* now E[s] == mean */
  293.         E[b] -= d;
  294.         if (E[b] < mean) {
  295.             push_stack(Smalls,b); /* no longer big, join ranks of the small */
  296.         }
  297.         else if (E[b] > mean) {
  298.             push_stack(Bigs,b); /* still big, put it back where you found it */
  299.         }
  300.         else {
  301.             /* E[b]==mean implies it is finished too */
  302.             (g->A)[b]=b;
  303.             (g->F)[b]=1.0;
  304.         }
  305.     }
  306.     while (size_stack(Bigs) > 0) {
  307.         b = pop_stack(Bigs);
  308.         (g->A)[b]=b;
  309.         (g->F)[b]=1.0;
  310.     }
  311.     /* Stacks have been emptied, and A and F have been filled */
  312.  
  313.     
  314. #if 0
  315.     /* if 1, then artificially set all F[k]'s to unity.  This will
  316.      * give wrong answers, but you'll get them faster.  But, not
  317.      * that much faster (I get maybe 20%); that's an upper bound
  318.      * on what the optimal preprocessing would give.
  319.      */
  320.     for (k=0; k<Kevents; ++k) {
  321.         (g->F)[k] = 1.0;
  322.     }
  323. #endif
  324.  
  325. #if KNUTH_CONVENTION
  326.     /* For convenience, set F'[k]=(k+F[k])/K */
  327.     /* This saves some arithmetic in gsl_ran_discrete(); I find that
  328.      * it doesn't actually make much difference.
  329.      */
  330.     for (k=0; k<Kevents; ++k) {
  331.         (g->F)[k] += k;
  332.         (g->F)[k] /= Kevents;
  333.     }
  334. #endif    
  335.  
  336.     free_stack(Bigs);
  337.     free_stack(Smalls);
  338.     free((char *)E);
  339.  
  340.     return g;
  341. }
  342.  
  343. size_t
  344. gsl_ran_discrete(const gsl_rng *r, const gsl_ran_discrete_t *g)
  345. {
  346.     size_t c=0;
  347.     double u,f;
  348.     u = gsl_rng_uniform(r);
  349. #if KNUTH_CONVENTION
  350.     c = (u*(g->K));
  351. #else
  352.     u *= g->K;
  353.     c = u;
  354.     u -= c;
  355. #endif
  356.     f = (g->F)[c];
  357.     /* fprintf(stderr,"c,f,u: %d %.4f %f\n",c,f,u); */
  358.     if (f == 1.0) return c;
  359.  
  360.     if (u < f) {
  361.         return c;
  362.     }
  363.     else {
  364.         return (g->A)[c];
  365.     }
  366. }
  367.  
  368. void gsl_ran_discrete_free(gsl_ran_discrete_t *g)
  369. {
  370.     free((char *)(g->A));
  371.     free((char *)(g->F));
  372.     free((char *)g);
  373. }
  374.  
  375. double
  376. gsl_ran_discrete_pdf(size_t k, const gsl_ran_discrete_t *g)
  377. {
  378.     size_t i,K;
  379.     double f,p=0;
  380.     K= g->K;
  381.     if (k>K) return 0;
  382.     for (i=0; i<K; ++i) {
  383.         f = (g->F)[i];
  384. #if KNUTH_CONVENTION
  385.         f = K*f-i;
  386. #endif        
  387.         if (i==k) {
  388.             p += f;
  389.         } else if (k == (g->A)[i]) {
  390.             p += 1.0 - f;
  391.         }
  392.     }
  393.     return p/K;
  394. }
  395.